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In principle, the generalized master equation can be used to efficiently compute the macroscopic 
first passage time (FPT) distribution of a complex stochastic system from short-term microscopic 
simulation data. However, computing its transition function matrix, T(r), from such data can 
be practically difficult or impossible. We solve this problem by showing that the FPT moment 
generating function is a simple function of the (easily computable) Laplace transform of the local 
FPT distribution matrix. Physical insight into this relationship is obtained by analyzing the process 
of steady-state relaxation. 
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Computing the macroscopic transition rates and first 
passage time (FPT) distributions of complex systems 
(e.g., proteins) by computational simulations using mi- 
croscopic equations of motion can be expensive, if not im- 
possible, particularly when the systems are large and/or 
when transitions are rare. In principle, computation 
can be simplified and comprehension can be enhanced 
by coarse-graining the microscopic equations to a macro- 
scopic generalized master equation [l| 

^ = S(t)P(0) - £ T(r) ■ P(t r)dr , (1) 

with P{— oo) = 0. This describes the time-evolution of 
P{t), the ensemble occupation number N- vector defined 
over states s, each corresponding to a subregion of the mi- 
croscopic phase space, after injection of systems at time 
t = 0. r(r) is the N x N matrix of transition func- 
tions, which includes memory effects. Corresponding to 
conservation of probability and to causality, 

l-r(r) = 0, rV(r)<0 (s^s'), (2) 

(where 1 is the A^-vector with all components equal to 
I)- 

The task is to determine the rate, or more generally 
the FPT distribution, of transitions from an initial state 
i to a final state /. So as to examine first passage times, 
we make / an absorptive state: 

r(r) ■ i f = , (3) 

where i s denotes the basis vector which has component 
s equal to 1 and all other components 0. Then the FPT 
distribution is 

ip(i)=dPf{i)/dt forP(0)=€ i5 (4) 

the mean FPT (MFPT) is the first moment {(rip)), where 

/>oo 

<(T fe ^»= / rVWrfr, 
Jo 



and the transition rate is ({rip))^ 1 . We assume that / is 
the only absorptive state and that the system is ergodic, 
so P(oo) — if. Therefore, 

/>OC 

/ V {T)dT = Pf(00) = 1. (5) 

Jo 

Defining the Laplace transform of g{r) as g(u) = 
J a e~ UT g{r) dr , the FPT moment generating function 
is 

L ~^r— = *»(-«) • 

fc=0 

We assume (as is true in most cases of interest) that f(r) 
decays faster than e _amaxT as t — > oo for some positive 
owx, so that if (—a) is analytic in a neighborhood about 
and can be differentiated to yield all moments. (This 
assumption is not essential, but simplifies the discussion. 
If it is not true, the discussion will still be valid for the 
finite moments.) 

In some cases, T(r) can be defined from first principles. 
However, for complicated systems such as proteins it 
must be computed from microscopic simulations. In such 
cases relatively short (compared to the MFPT) molecular 
or Langevin dynamics simulations can be used to deter- 
mine —K S)S i(t) (s =/= s'), the local FPT distribution. This 
is the probability density that, after arriving at state s', 
a system waits for an interval t before first leaving and 
that it goes to s. The diagonal elements of the N x N 
matrix K(t) are defined by K s , s (t) = - J2 S ^ S ' k s,s>(t), 
so like T(t), K(t) satisfies 

1-X(r)=0, X s>s ,(r)<0 K(T)-e f = 0, 

and also, by its definition and the assumption that / is 
the only absorptive state, satisfies 

poo 

/ K s , s ( T )dT=\ (s^f). 

Jo 
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The Laplace transforms of T(t) and K(t) are related by 



r(u) = uK{u) ■ {I - BiaglK(u)]}- 1 
K(u) = r^-juJ + Diaglr^)]}- 1 , 



(6a) 
(6b) 



where / is the identity matrix and Diag[A] denotes the 
matrix of diagonal elements of A Q[. However, even 
though Eq. I|6a|l can be used to determine T(u) from 
K{u), the inverse Laplace transform required to deter- 
mine r(r) can be difficult, if not impossible, to compute. 
This precludes the computation of <p(t) by integration of 
Eq. 0) when only K(t) is known. 

Faradjian and Elber 0] have recently provided a so- 
lution to this problem by showing that Eq. can be 
reformulated using either a coupled set of QK equations 
or a single set of (pseudo-)Markovian PJ equations that 
allow dP(t)/dt to be integrated in t using only K (r), not 
r(r). However, this procedure is inefficient since it re- 
quires modelling the full functional dependence of K(t), 
even though most of the relevant information is contained 
within only a few of its lower moments. 

Here we discuss how to compute the FPT moments di- 
rectly from the moments of T(t) or K(t) without need for 
an alternative formulation and at reduced computational 
cost. In addition, we show how physical insight into the 
relationships between the moments emerges from analy- 
sis of the physical process of steady-state relaxation. 

Results — The Laplace transform of Eq. 0) with 
P(0) = e z is 



uP(u) 



r(o) • p(u) 



with solution P(u) = [ul + r(ii)] _1 • e$, where I is the 
identity matrix. Combining this with Eqs. and Qfiajl 
gives expressions for the FPT moment generating func- 
tion in terms of either T(u) or K{u): 

(p(—a) = — aif ■ [r(— a) — al}^ 1 ■ e< (7a) 
= if ■ [I + fi-a)}- 1 -ii , (7b) 

where A denotes the matrix of off-diagonal elements of A 
0- [The singularity in Eq. (|7a|l at a = is removable.] 

To gain insight into the physical significance of these 
relationships, we replace Eq. with the equation for 
steady-state relaxation — the situation in which systems 
are continuously injected into the initial state at an ex- 
ponentially decreasing rate exp(— at) beginning at t = 0. 
Eq. JIJ becomes 



dP(t) 
dt 



/ T(r) • P(t - r)dr , 
Jo 



(8) 



where 8(t) is the Heaviside step function. If a = 0, as 
t — > oo Eq. (JSJ) describes steady-state flow, and we ex- 
pect that the occupation numbers in all the intermediate 



states, P s (t) (s 7^ /), will approach a constant steady- 
state solution. When a > we expect that the interme- 
diate state occupation numbers will will decay exponen- 
tially as t — > oo. The key point is that the asymptotic 
steady-state solution is easy to compute and determines 
the FPT generating function exactly. 

Since Eqs. (JIJ and © are linear, the systems injected 
at different times will transition to / independently, so 
dPf /dt during steady-state relaxation is obtained simply 
by integrating Eq. Q over the incoming flux: 



dPfjt) 
dt 



-a(i-T)^ T ) dr _ 



(9) 



That is, <p{r) determines the fraction of the flux intro- 
duced at time t — r, e~ a ^~ T \ that arrives at / at time 
t. As long as t is larger than the support of (p(r), we can 
extend the upper limit of the integral to oo, and Eq. @ 
can be rewritten as: 



(p{-a) = / <p{r)e aT dT 



, dP S (t) 
dt 



(t»0) 



(10) 



Thus, solving Eq. (JSJ for P{t) when t 3> determines 
the FPT generating function. 

We must take care because T(r) is singular: Eqs. 
and J2J imply that it has left null-vector 1 and corre- 
sponding right null- vector if. Thus, it is convenient to 
decompose Eq. (JHJ by projecting it into the null-space 
and its bi-orthogonal space using the idempotent, asym- 
metric projection operators V and Q: 



P = e/®1, Q = I-V 



and expanding 



P(t) = i f P t ot(t)+P(t) , 



(11) 



where Ptot(t) = 1-P(t) and P(t) = Q-P(t). Substituting 
this into Eq. JHJ yields the two independent equations 



dPtotjt) 
dt 

dP(t) 
dt 



= e- at 9(t) 

= e-^ewiii-if) 



T{t) ■ Pit - t) dr, 



The first equation [along with boundary condition 
P(— oo) = 0] implies that 



Ptat(t) 



1 - e 



-6{t) 



(12) 



We solve the second equation in the asymptotic regime by 
testing the guess that P(t) = P exp(-at) (t 3> 0), where 
1 • P = 0. Substituting this in, factoring out exp(— at), 
and taking the limit t — > oo gives 

-aP = (€i-if)-f(-a)-P (13a) 
=> P = Q- [f (-a) - al}- 1 ■ - i f ) . (13b) 
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Combining Eqns. I|10H13[1 with Eq. I|6a|l gives the gen- 
erating function as 

<p{-a) = 1 - ae f ■ Q • [f {-a) - al}- 1 ■ Q ■ l % , (14) 

which is equivalent to Eqs. J7J). In this form it is evident 
that the zeroth moment is always </5(0) = 1, consistent 
with Eq. (J3J). The higher moments are obtained by ex- 
panding in powers of a about a = with care for the 
degeneracy of T(0). For example, the first moment in 
terms of F is 



d(p{— a) 



da 



(1 



r(o) s 



(15) 



where A 3 denotes the spectral (or Drazin) pseudoinverse 
0. 

Eq. I|15(l has a simple interpretation: The steady-state 
(a = 0) flux must be equal to the total number of sys- 
tems in transit divided by the MFPT. Thus, when the 
flux equals 1 (as in this calculation), the MFPT must 
be equal to the sum of the occupation numbers in all the 
intermediate states, Y^sjtf ^s(*)- Eqs. (fTTfl and (|13bl) im- 
ply that the rhs of Eq. (|15fl equals (1 — if) ■ P, which is 
identical to this sum. We also note that Eq. Q15[l is the 
same result that would be obtained if we ignore all mem- 
ory effects and approximate T(r) ss S(r) J Q r(r)dr; this 
is equivalent to replacing the generalized master equa- 
tion with a regular master equation having T — J T(r). 
Differences between the moments of these two equations 
only appear in higher order. 

When r(it) has to be determined from K(u), we can 
combine Eqs. lj7a|l and (|T3|l [or directly Eq. i|7b|l ] to get 



((rip)) = l-Dmg(((TK)))-[I + K(0)}- 
= 1 • Diag(«r^))) ■ K(0) s ■ it . 



(16) 



Thus, the MFPT is completely determined by the zeroth 
moments of K (r) and the first moments of its diagonal 
elements, which are the mean incubation times of the 
individual states. Only these TV + Nj^ moments [where 
Njf is the number of non-zero off-diagonal elements in 
K(t)] must be determined from the molecular simulation 
data to compute the MFPT. 

Not only is Eq. lfTo|l simple to compute, but we expect 
that it will require less simulation data than direct in- 
tegration to compute the MFPT to the same statistical 
accuracy: Integrating dP(t) jdt requires determination of 
a complete discretized representation of K(t): i.e., the 
(N+N^xTr/h values oiK{nh) s , s > for n = 0, 1, . . . T T /h, 



where h is the smallest relevant time scale and T T is the 
support of K(t). Each value must be determined by 
counting the number of simulation events that make the 
first passage from s' to s in time nh < r < (n + l)h. 
Since each number is a small fraction of the total number 
of simulation events, the statistical error of its determi- 
nation will be large relative to that of the zeroth and first 
moments. Therefore, direct integration will require more 
molecular simulation data than Eq. Ijl6(l to achieve the 
same statistical error. 

Analogous expressions for the higher FPT moments 
can be determined by computing additional values of 
tp(— a) near a = and numerically differentiating or 
by analytically expanding Eq. Q) or l|14|) in terms of 
the ((T k K)) to obtain expressions analogous to Eq. (|lfc>|> 
The latter procedure can be used to compute the fi- 
nite moments even when higher moments are infinite and 
(p{— a) is not analytic at a = 0. 
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